!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate point on exit Mach line for a supersonic  
! wind tunnel nozzle dedign with Mexit=Md by MOC.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-07-15
! revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : Exitline
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine Exitline(alpd,stepd)
!          4()
!      2
   use VariableDef
   implicit none
!   given variables
!   real*8::x2,y2,Vu2,Vv2
   real*8::alpd,stepd
   integer::icor,i
!   temp variables
   real*8::Vup,Vvp,yp,Vp,thetap,ap,Mp,lp,Qp,Rp,Sp,Tp
   real*8::x4f,Vu4f
!   calculate variables
!   real*8::x4,Vu4,V4,p4,T4,rou4,M4
   icor=6
   Vu4f=0.0
   y2=abs(y2)
   Vv1=abs(Vv1)
   Vup=Vu1
   Vvp=Vv1
   yp=y1

   x4=x2+stepd*cos(alpd) !/tan(alpd)
   y4=y2+stepd*sin(alpd)
   Vv4=0.0   ! streamline parallel to x-axis

   do i=1,icor
      Vp=sqrt(Vup**2+Vvp**2)   ! Vup
      thetap=atan(Vvp/Vup)   ! 0.0
      call Thermo(Vp,p4,rou4,T4,ap,Mp)
      
      lp=tan(thetap+asin(1/Mp))   ! tan(asin(1/Mp))
      Qp=Vup**2-ap**2
      Rp=2*Vup*Vvp-Qp*lp   ! -Qp*lp 
      Sp=delta*Vvp*ap**2/yp   ! 0.0
      
      Tp=Sp*(x4-x2)+Qp*Vu2+Rp*Vv2
      
      Vu4=Tp/Qp

      if (abs(Vu4-Vu4f) <= ev) then
         exit
      end if

      Vup=(Vu2+Vu4)/2
      Vvp=(Vv2+Vv4)/2
      yp=(y2+y4)/2
      
      Vu4f=Vu4
   end do
   V4=sqrt(Vu4**2+Vv4**2)
   call Thermo(V4,p4,rou4,T4,ap,M4)
   
   if (wallflag==-1) then
      y1=-y1
   end if
   
end subroutine Exitline
